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Due to a limited number of accelerometers available for use, the shock 
trial for the DDG-51 class destroyer provided a spatially incomplete set of time 
history data. However, a visualization of the shock response of the entire ship 
is desired. To this end, finite element model reduction methods are employed 
to provide a transformation matrix which is used to expand this relatively 
small collection of data into the same number of degrees of freedom as the 
finite element model. Using this expanded set of time histories, it is possible 
to animate the transient response of the structure as a whole. 

This approach is investigated using computer-simulated transient 
response data from a finite element model of a flat plate. The use of static and 
dynamic reduction methods are explored in the creation of the 
transformation matrices required for the visualization of the expanded data. 
The animations are assessed based on a quantitative comparison with the 
full-order transient model response. 
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I. INTRODUCTION 



To gain information concerning the transient response of a complex 
structure under an arbitrary loading, an analysis of vibration response time 
history data is required. Unfortunately, a continuous system has an infinite 
number of degrees of freedom from which data can be extracted. Taking 
measurements at all of possible locations on the structure is obviously not 
feasible due to constraints on time, resources, and money. Because of these 
constraints, accelerometers are placed at a few specific locations on a structure 
and the data is then analyzed. While this information is quite useful, the 
collection of data still does not provide a clear picture of the response at the 
many locations where data was not taken. More importantly, since each piece 
of data is for a specific location, it takes quite a bit of imagination to get an 
understanding of how the structure as a whole is responding to the loading. 

Due to the limits on the number of accelerometers available, the shock 
trial conducted on the DDG-51 class destroyer provided a spatially incomplete 
set of test data. In addition to the information collected during the test, 
knowledge of the transient response at locations on the ship where data was 
not taken was desired. By expanding this collection of test data into a spatially 
complete set, the transient response could be visualized in a computer 
simulation. 

This thesis investigates the strengths and weaknesses of expanding 
spatially incomplete test data using finite element model reduction methods. 
These reduction methods are used to create a transformation matrix which 
makes it possible to expand a relatively small collection of simulated test data 
into a set that corresponds to a much larger number of degrees of freedom. 
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The transformation matrix is created from a finite element model of the 
structure whose system mass and stiffness matrices have been partitioned in 
accordance with the locations of the accelerometers on the actual structure. 
The final collection of data becomes the time histories of all the coordinates 
in the full-order finite element model and has the same number of degrees of 
freedom. 

The full-order set of data is then animated to create a real-time 
visualization of the dynamic response of the structure as a continuous 
system. Every column in the expanded matrix is a representation of the 
dynamics of each node in the finite element model at a specific time interval. 
With this information, a very clear understanding of the vibration response 
of a complex structure is possible. 

There are two different reduction methods used in this thesis to create 
the transformation matrices. The first is the static, or Guyan, reduction and 
the second is the Improved Reduced System (IRS) reduction. The IRS 
method, while initially more computationally expensive, provides an 
improvement because it approximates the inertia forces that are present in a 
dynamic problem. A sensitivity analysis for each simulation is provided to 
highlight the accuracy of both the static and the dynamic reduction methods. 
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II. THEORY 



A. MATRIX PARTITIONING 

In the physical coordinate system, the expansion of the collected time 
history data can be expressed in the following way: 




(1) 



where {x a } is a vector containing the dynamic response data collected from 
the structure, {x a } is the response at coordinates other than the data collection 

points, and [T] is a non-square transformation matrix used in the expansion 
of {x a } to the full-order response set, and [I] is the identity matrix. 

The subscript 'a' refers to the analysis set of model coordinates, which 
from here on will be called simply the 7 a-set 7 . This is the set of coordinates in 
the model that corresponds to the locations on the structure where data was 
taken. Conversely, the subscript 'o' refers to the omitted set. This '0-set 7 refers 
to all subsequent model coordinates where data is not taken. By operating on 
the 'a-set 7 time history data set with the proper transformation matrix, the 7 o- 
set 7 dynamic response is calculated. 

Generally, due to the size of complex models, the a-set system will be 
much smaller than the o-set system. The model of a complex structure will 
normally contain a huge number of degrees of freedom while the number of 
accelerometers is relatively small in comparison. 

B. STATIC TRANSFORMATION MATRIX 

By partitioning the stiffness relation, f = kx, the static transformation 
matrix can easily be found. The equation becomes: 
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( 2 ) 



■K« K Jfxl ffl 

K« K^JjxJ jfj 

By taking care to partition the equation such that there are no excitations 
applied at the o-set coordinates, then {f 0 } can be set equal to {0}. Now the 

relation between the o-set and a-set response can be expressed as: 



kMxXjk} 



The static transformation matrix is then [Ref. 1], 



T = 

star 






(3) 



(4) 



C IRS TRANSFORMATION MATRIX 

The derivation of the IRS reduction transformation matrix begins with 
the partitioned equation of motion for a vibrating system [Ref. 2]: 



M 

M 



aa 

oa 




K u 

K„ 




(5) 



Assuming a harmonic motion of frequency, Q, the acceleration term in 
equation (5) can be replaced with x = -Q 2 x. Again partitioning such that there 
are no excitations at the o-set coordinates, we arrive at the exact structural 
dynamics reduction relation [Ref. 3]: 
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( 6 ) 



M = [I - a ! Ol..]"[-K;X, + n ! K^M„]{x.} 

(a) (b) 

Ideally, the desired transformation matrix will be frequency independent and 
Eq. (6) is dependent on the driving frequency by the £2 2 term. Truncating the 
binomial expansion of Eq. (6a) after the terms that include Q. 2 yields: 

K}[-k;;k„ + n 2 dCM„„ - k^mjOc. )]{*.} (7) 

The frequency dependency is then eliminated by the following approximation 
of the acceleration: 



n 2 {x.} = [M, ra ]-'[K,„] (8) 

where the statically reduced mass and stiffness matrices are given by the 
following [Ref. 2,3]: 



M„,=Ti,MT„ (9a) 

K„„=TI»KT„, (9b) 

Substitution of Eq. (8) into Eq. (7) provides: 

{x„}[-K^K m+ [K-jM,.-K^M„K^K„][M„„]-'[K„„]]{ X ,} (10) 

which is the IRS reduction relationship between the a-set and o-set 
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coordinates. The IRS reduction transformation matrix can then be written as 
[Ref. 2,3]: 



T = 

1 1RS 



-K^K_ +[K^M„ 



I 



(ii) 



D. VISUALIZATION MATRIX 

After operating on {x a } with the transformation matrix, the response 
for the entire structure is partitioned in the form: 



{x(r)} = 



f*.(0 1 


K(0} • 


• Kco} 


k(oj 


.Ka,)} • 


• K(0}_ 



( 12 ) 



where the columns in the matrix are a representation of the dynamic 
response at each node in the model at a specific instant in time. However, this 
partitioning does not correspond to the nodes in the model, therefore Eq. (12) 
must again be partitioned such that it is in line with the nodes in the model. 
To get a real-time visualization of the structure, simply march through the 
matrix taking snapshots of the model at each time interval. 
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III. FINITE ELEMENT MODEL FORMULATION 



The finite element method used in the creation of the plate model for 
this thesis is based on the shear deformable displacement formulation. The 
element chosen is the four-noded quadrilateral plate element. Each node in 
the model has three degrees of freedom which are the transverse 
displacement, w , and the two rotations about the midplane, G x and 9 y . Using 

the following bilinear isoparametric shape functions [Ref. 4], 



^1(^77) = ^-(1 - ^)(1 - 77) 


(13a) 


^O+QO- 7 *) 


(13b) 


H A&'n) = \( l + Z )( 1 + 7 l) 


(13c) 


90 + 17 ) 


(13d) 


the transverse displacement and slopes are interpolated as: 




n 

w= X // .(^ r iK 

/=1 


(14a) 


;=i 


(14b) 


e r =im.n)(e,l 

i —\ 


(14c) 


These isoparametric shape functions are defined in terms of a 
domain -1 < £ < 1 and -1 < T] < 1 . 


normalized 



A. INTEGRATION TECHNIQUE 

As previously stated, the shape functions used in this model are for the 



bilinear isoparametric element. This type of element was chosen in order to 
make the finite element model general in the sense that it can be applied to 
any plate geometry. The elements in the physical coordinate system are 
mapped into the isoparametric coordinate system by the following relations: 
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x =^Z H i(^) x i 

1=1 


(15a) 


y = 'L H i{Z’'n)y i 


(15b) 



1=1 



where x and y are physical coordinates corresponding to the nodes in the 
element. Gauss-Legendre quadrature is used to perform all integrations. 

B. ELEMENT MATRIX 

1. Elemental Stiffness Matrix 

Without providing the derivation, the elemental stiffness matrix, [K e ], 
for plate bending is [Ref. 4]: 



[K e ] = ~ J + kh j BjDjBjtlQ 









(16) 



where h is the thickness of the plate, k equals — and is the shear energy 
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correction factor. Of is the two dimensional element domain, and 



[B*] = 



dH, 

dx 


0 

dH, 


0 


dH 2 

dx 


0 0 
dH, n 


dH, 

dx 


0 


1 

dy 


0 


0 


— 2- 0 

dy 


0 


dH, 


dH, 


0 


dH 2 


d Jk 0 


dH, 


dy 


dx 


dy 


dx 


dy 



0 

dH, 

dy 

dH, 

dx 



0 *1 

dx 

0 0 
dH 



0 0 



dH, 

dy 



0 



0 



« 3H < 0 



dy dx 



(17) 
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0 



0 



0 



[B s ] = 



-Hi 

0 



-H x 



<?y 



-h. 



0 -H, 



<?h 2 

* 

dH 2 



-Hi 



0 



-H, 



<^3 

<2r 

£^3 

<?y 



-H. 



0 

-H, 



<9H 4 

dx 
dH 4 

(18) 



The constitutive equation for the plane stress condition is: 




1 v 
v 1 
0 0 



0 

0 

1 — v 
2 



while the constitutive equation for shear is. 




0 

G 



(19) 



(20) 



G is the shear modulus, E is the modulus of elasticity, and v is Poisson's 
ratio. One thing of importance to be noted about Eq. 15 is that the shear energy 
becomes dominant over the bending energy as the thickness of the plate 
becomes small relative to the side length. This is called shear-locking. To 
account for this problem, a reduced integration technique is used. The 
bending term is integrated exactly using 2x2 Gauss-Legendre quadrature while 
the shear term is under-integrated using lxl Gauss-Legendre quadrature. 

2. Elemental Mass Matrix 

The elemental mass matrix is found by simply integrating the properly 
weighted shape functions. The general equation for determining this matrix 
is: 



9 



( 21 ) 




In Eq. (19), P is the density of the plate element, A is the elemental area, and 
[ZZ] is the same shape functions as in Eq. (12). 2x2 Gauss-Legendre quadrature 
is used to perform all integrations for the element mass matrix. 
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IV. COMPUTER SIMULATIONS 



The investigation of expanding time history data was conducted with 
the aid of computer simulations. The dynamic response at each node of a 
plate finite element model was used to simulate the 'actual' solution. From 
this baseline model, an a-set system of time histories was chosen which is 
assumed to be the data taken from an actual structure. The system matrices in 
the finite element model were then partitioned in order to create the 
transformation matrix required to expand the a-set time histories to full- 
order. By comparing the dynamic response of the finite element model to the 
response found by expanding the a-set system, the validity of the process 
could be studied and verified. 

A. EXAMPLE (1): TWO DOF IN THE A-SET 

The initial simulations are conducted on a square, flat plate (see Figure 
1) comprised of twenty-five elements and thirty-six nodes. Degrees of freedom 
42 and 87 are assumed to be the 'actual' response and therefore make up the a- 
set system. Conversely, the o-set system is comprised of all other DOF. The 
two DOF in the a-set system correspond to the transverse motion at nodes 14 
and 29. The external force is a blast function (see Figure 2) applied at node 
twenty-one. Each external spring has a stiffness equal to 1000 lb-in and is 
attached to the ground. 

Two separate simulations are performed using this plate and a-set 
system. The first simulation uses the static transformation matrix to expand 
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the a-set and the second uses the IRS. Both are then compared to the full- 
order plate response. Unfortunately, it is impossible to show the visualization 
on paper so the comparison will be on a per DOF basis. 




Figure 1 Plate Model 
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Figure 2 Forcing Function 
1. Full-Order Plate Results 

The dynamic response of the full-order plate model was solved by 
transforming the general equation of motion, 

[Mp} + [^]{x} = {/(0} (22) 

into modal coordinates by the following relation: 
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Mt)} = [<t>]{q(t)} 



(23) 



in which (q(t)} is the modal displacements (as a function of time) and [<))] is a 
matrix of mass normalized modeshapes found from solving the typical 
structural dynamics eigenvalue problem: 



The modal displacements are found by numerically solving a convolution 
integral [Ref. 5]. The damping in the plate is assumed to be 2% and is inserted 
into the integral calculation. Once [q (t)} is found, Eq. (21) is then used to 
convert back to the physical coordinate system. 

Now, the time history at every node in the plate is known. These results will 
be used to compare the accuracy of the static and dynamic expansions of the a- 
set. 

2. Static Transformation Results 

Figures III, IV, and V are per-node samples of the full-order plate 
model response in comparison to the response found by the static expansion 
of the 2 DOF a-set. It is clear that the static transformation will provide an 
adequate approximation of the 'true' response at certain nodes. 

Unfortunately, Figure 3 proves that the difference between the 'true' solution 
and the expanded set can be greater than 50%. 




(24) 
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Figure 3 Static Response at Node 3 
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Figure 4 Static Response at Node 17 
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Time History at Node 27 
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Figure 5 Static Response at Node 27 

Having errors of this magnitude can lead to very misleading results in 
the final visualization. It will provide an idea of how the structure vibrates 
but the animation is distorted and skewed. If a more accurate representation 
is the requirement, then other steps must be taken. Obviously, one way to get 
a more accurate solution is to increase the size of the a-set system. An 
example solution of the larger a-set is provided later. The drawback to this is 
more data must be taken which is oftentimes impossible due to limits on the 
number of accelerometers available. Another way to change the accuracy of 
the solution is to alter the locations where the data is actually taken on the 
structure. An example is also provided later in this chapter of greatly 



/■ *\ 

/ \ 



/ \ 

/ > 



_ Static: - 
_ Model: -. 



\ / 
\ / 



16 



decreasing the accuracy of the static solution by choosing an a-set system 
badly. 

a. Error Analysis 

To understand the differences in the full-order and static plate 
solutions without providing the time histories at every node, a sensitivity 
analysis was conducted. 



Sensitivity Analysis Between Static and Full-Order Response 




Figure 6 Error of Static Reduction Method 
note: Figure is orientated such that the 0/0 point is node 1 



The largest difference between the full-order plate response and 

the expanded time history set was found and plotted. This analysis is 

essentially a 'worst-case' scenario because the code searched through time to 

17 



find the absolute maximum error at each node and these errors did not 
always occur at the same instant in time. 

The maximum differential appears to be approximately 0.1 
inches, but the movement of the plate is only about 0.3 inches. Under these 
circumstances, an error of this magnitude is quite significant, and therefore 
the static reduction probably should not be used. 

3. IRS Transformation Results 

Using the existing plate model and forcing function (see Figures 1 and 
2), the solution for the IRS expanded data is provided as before. The same 
nodes are given as samples here: 
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Figure 7 IRS Response at Node 3 
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Figure 8 IRS Response at Node 17 



Time History at Node 27 




Figure 9 IRS Response at Node 27 
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These figures show that the IRS expanded data is much more accurate 
in comparison to the static transformation. Even with an a-set of only 2 DOF 
out of the possible 108, the IRS transformation provides a very good solution. 
Clearly, the correction in the IRS transformation matrix due to inertial forces 
will improve the accuracy of the solution immensely, 
a. Error analysis 

A sensitivity analysis for the IRS expanded data is also provided: 



Sensitivity Analysis Between IRS and Full-Order Response 




Figure 10 Error of IRS Reduction Method 
note: Figure is orientated such that the 0/0 point is node 1 

The maximum difference in Figure 10 is less than 0.02 inches. 
This accuracy is an order of magnitude better than the static reduction. The 
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approximation of the inertial forces inherent in the IRS transformation 
matrix appears to work quite well. 

B. EXAMPLE (2): FOUR DOF IN THE A-SET 

This example shows that the accuracy in both solutions can be greatly 
improved by increasing the size of the a-set. By merely increasing the size of 
the a-set by two, it will be seen that the improvement in accuracy is at least an 
order of magnitude. As mentioned earlier, although desirable, increasing the 
number of data points may not be feasible. 

For this simulation, the forcing function used is the same as in 
Figure 2. The plate model, shown in Figure 1, is changed such that the data is 
instead extracted at nodes 8, 11, 26, and 29. In the transverse direction, these 
nodes correspond to degree of freedom 24, 33, 78, and 87. All other variables 
in the simulation are unchanged. 

On a per node basis, the results are very hard to comprehend because 
the graphs are almost exactly on top of each other. For the IRS expanded data, 
no difference can be detected. Because of this, only the sensitivity analysis is 
provided. 

1. IRS and Static Transformation Error Analysis 

By viewing these graphs, it becomes clear how the accuracy can by 
greatly increased by adjusting the number of DOF in the a-set system. The 
greatest difference in inches of the response for the static solution is 0.01. 
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Response Differential (in.) 



Sensitivity Analysis Between Static and Full-Order Response 




Figure 11 Error of Static Reduction Method 
Sensitivity Analysis Between IRS and Full-Order Response 




Figure 12 Error of IRS Reduction Method 
note: Figures are orientated such that the 0/0 point is node 1 
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The difference for the IRS solution is much smaller at 0.00015 in. As expected, 
the IRS reduction method remains the more accurate of the two, but both 
methods provide very good solutions. If the visualization could be shown 
here, one would not be able to discern a difference between the two 
animations and the 'true' solution. 

C EXAMPLE (3): A-SET INCLUDES CONSTRAINED NODES 

For this simulation, the accelerometers are assumed to be located at 
nodes which are constrained to ground through external springs. For the 
plate, these nodes are numbered 1, 6, 31, and 36. The corresponding DOF's in 
the a-set system are 3, 18, 93, and 108, respectively. 

1. Static Transformation Results 

Figures 13 and 14 are a sampling of the results from this particular 
simulation and they show a very interesting situation to consider. These time 
history samples clearly indicate that the static transformation matrix will not 
provide an accurate solution for this particular a-set. On closer inspection of 
the two previous graphs, it is seen that the static transformation provided the 
exact same solution at nodes 10 and 29. In fact, although not shown, this is the 
exact response calculated at every node in the plate, and, when placed in 
animation, the plate exhibits rigid body motion in the transverse direction. 

The reason for this rigid body motion can be attributed to the location 
of the a-set system. By taking time history data at only those nodes 
constrained to ground through springs, the transformation matrix imparts 
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the strain energy of the springs to the system and ignores any of the plate's 
internal strain energy. In other words, the response calculated in the o-set is 
effectively 'zeroed out' by the transformation matrix. By discarding the 
motion of the individual o-set nodes, the plate moves with the springs as a 
rigid body. A detailed explanation of this behavior is provided in the next 
chapter. 

a. Error analysis 

The inaccuracy of the solution becomes clear with the sensitivity 

analysis. 



Sensitivity Analysis Between Static and Full-Order Response 




Figure 15 Error of Static Reduction Method 
note: Figure is orientated such that the 0/0 point is node 1 



25 



Obviously, rigid body motion is an inaccurate solution to this 
situation. As expected, the nodes at the four comers exhibit no error while the 
center of the plate has the greatest. From Figure 13, the maximum amplitude 
of the plate motion is about 0.3 in. which is the same as the maximum error 
shown here. Example (2) proved that a very accurate solution can be found 
using the static reduction method when the a-set contains four DOF. This 
example shows that care must be taken when using the static reduction 
method that the accelerometers are not placed at constrained nodes. 

2. IRS Transformation Results 

The problems using the static reduction method are not exhibited by 
the IRS method. Due to the correction in the derivation of the 
transformation matrix due to the inertial forces, the motion of the plate is 
taken into account. For this reason, the actual response of the plate is 
calculated. 

The following time history samples prove the usefulness of the IRS 
reduction method when data must be taken at nodes constrained to ground 
through a spring: 
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Figure 16 IRS Response at Node 10 



Time History at Node 29 




Figure 17 IRS Response at Node 29 
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a. Error analysis 

Figure 18 has the same basic shape as Figure 15, but the 
maximum error this time is only 0.02 in. The inertial force correction to the 
IRS transformation matrix enables a very accurate solution to be found even 
though the strain energy of the plate is ignored. For this reason, using the IRS 
method allows for fewer limitations in the location of the a-set system. 



Sensitivity Analysis Between IRS and Full-Order Response 




Figure 18 Error of IRS Reduction Method 
note: Figures are orientated such that the 0/0 point is node 1 
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V. LESSONS LEARNED 



The three examples provided in the previous chapter are only a few of 
the simulations conducted in the investigation of expanding spatially 
incomplete time histories to full-order. Several other situations requiring 
investigation were checked and will be addressed here. This chapter 
highlights the lesson's learned from all the computer simulations without 
providing the detail seen previously. 

A. ACCURACY DEPENDS UPON THE NUMBER OF DOF IN THE A-SET 

As expected, the number of DOF in the a-set system greatly affects the 
accuracy of the solution. The increase in accuracy seen in changing the a-set 
system from Example (1) to Example (2) is excellent proof of this. Actually, a 
dramatic improvement could be seen by increasing the number of DOF in the 
a-set to three. Interestingly, independent of the size, shape, or number of DOF 
in the model, a system of only one DOF in the a-set did not provide a viable 
solution for either method. Therefore the minimum number of DOF in the 
a-set is restricted to two. Simulations were not run for situations where the 
number is greater than four because, as seen in Example (2), an a-set 
containing four DOF expands into a solution that is essentially exact. 

B. LOCATION OF THE A-SET IS CRITICAL 

In addition to the number of DOF, the accuracy of the solution can be 
greatly affected by the location of the a-set. Two general cases have been found 
that should be taken into account. The first is to pick an a-set location that is 
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certain to impart strain energy into the transformation matrix in order to 
avoid the rigid-body modes. This restriction is of vital importance if the static 
reduction is chosen. The second is to ensure that the time history data is 
taken at points 'distant' from one another. 

Example (3) revealed a problem with the static reduction method when 
the nodes chosen for the a-set are constrained. As mentioned earlier, the 
reason is that the only strain energy imparted into the system is the stiffness 
of the springs while the internal strain energy of the plate is 'zeroed out'. 
Mathematically, this idea can be represented by first partitioning the stiffness 
matrix in the following manner: 
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where the superscript p is used to indicate the plate and s refers to the 
springs. Let x be a displacement vector that imparts no internal strain energy. 
This now defines the a-set and the o-set. Further partitioning of Eq. (23) yields: 
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Now, applying Eq. (24) to the general stiffness relation, f = kx, the following 
relation between the a-set and the o-set is derived much in the same manner 
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as Eq. (3): 



kk-KlKlM 



(27) 



By condensing the stiffness of the grounded springs out of the static 
transformation equation, the o-set exhibits the response of an unrestrained 
structure. The rigid body motion exhibited in Example (3) is a result of the 
stiffness matrix being partitioned in this manner. 

A problem that affects the accuracy of both the static and IRS reduction 
methods is choosing the a-set system at points that are 'too close' together. In 
this situation, common sense must be depended upon because the concept of 
'close' and 'far' depends on the size and geometry of the structure. Several 
simulations were conducted to study this and the best solution is simply to 
ensure the accelerometers are spread out to locations that best cover the entire 
structure. Having said this, it is important that the a-set include locations 
where the motion of the structure has the greatest amplitude. In the case of 
the plate model, the amplitude of the vibration is greatest at the center of the 
plate so the a-set should include time-history information from this area. A 
simulation was run where the a-set was comprised of edge nodes. The plate 
moves very little at its edges and subsequently the solution was inaccurate for 
both methods. Sound engineering judgment must be employed when 
choosing an a-set because it is critical to getting an accurate solution. 
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C ROTATIONS AS PART OF THE A-SET 



In general, using rotations, or a combination of rotations and 
translations as part of the a-set does not seem to alter the accuracy of the 
solution. The same rules apply for rotations as for the transverse 
displacements. A larger a-set still returns a more accurate solution and the 
location is still important. In this situation, however, being closer to the edge 
produced more accurate solutions. This makes sense because the edges are the 
locations of the largest rotations. Using rotations in the a-set is an advantage 
if applying the static reduction method because there is no longer a restriction 
on the use of nodes constrained through translational springs. 

While using rotations is not a problem computationally, care must be 
taken to ensure that the structure is rotating along the axis corresponding to 
the chosen DOF. If not, including these DOF in the a-set serves no purpose as 
no strain energy is imparted into the system and the plate will again exhibit 
rigid body motion. 

D. VISUALIZING A STRUCTURE 

Animation in MATLAB is very computer memory intensive. A 

simple structure such as the plate can tax the memory of most computers and 

cause the program to crash. As such, a very large structure, or one with many 

nodes, would be almost impossible to visualize. To aid in overcoming this 

problem, the code should direct the computer to reduce the size of the 

graphics window for the animation. Although this does not allow for as nice 

a visualization, memory is saved for computation instead of being allocated 
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to the graphics window. By doing this, the visualization can be of longer 
duration thereby resulting in a better understanding of the dynamics of the 
structure. 

A better tool to use in performing a dynamic visualization is I-DEAS. 
Unfortunately, certain versions will not allow the user to expand time 
histories. An outline can be found in the appendix that describes the steps 
needed to perform an animation in I-DEAS. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



The main objective of this thesis was to investigate the process of 
visualizing transient structural response by expanding spatially incomplete 
time history data to full-order. The response at every node of a full-order 
plate model was solved which provided a base model for comparison. From 
this collection of nodal responses, a few of the time histories were extracted 
which simulated the experimental data taken from an actual structure. After 
choosing the a-set, the finite element system matrices were partitioned to 
calculate the appropriate transformation matrix. The expanded data was then 
compared to the full-order response in the effort of evaluating the accuracy of 
the process. 

A. CONCLUSIONS 

The simulations have shown the following: 

1. To get an accurate* solution, the minimum number of DOF's allowed 
to make up the a-set system on a flat plate is two. Provided the 
accelerometers are properly placed, an a-set comprised of four DOF 
will produce a very accurate response for both the static and the IRS 
methods. 

2. The placement of the accelerometers is very important. If the static 

reduction method is chosen, the data must be taken at nodes certain 
to produce internal strain energy. Regardless of method, the data 
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must be extracted at locations 'far' enough from each other to get a 
sample of the entire structure, and the data collector should try to 
pick points on the structure where vibration is greatest. 

3. The same restrictions that apply to using translations in the a-set can 
be applied to rotations. 

4. The visualization process is very memory intensive therefore clever 

programming is required to maximize the memory available for 
graphics. 

Visualizing transient response provides a unique insight on the 
motion of a vibrating structure. A much better understanding of the 
dynamics is gained when a designer can visualize the motion of a spatially 
complete body instead of the time histories of a few select nodes. 

B. RECOMMENDATIONS 

This study outlines a few of the strengths and weaknesses inherent to 
the Guyan and IRS reduction methods in the attempt to animate transient 
response. While many simulations were explored, the investigation was 
rather limited in scope. In addition to the Guyan and IRS methods, there are 
many reduction methods that could be explored. To validate the process, 
experimental data must be loaded into the program and studied. A similar 
study should be directed towards visualizing the response of three 
dimensional structures. In addition, it would be useful to investigate 
animating transient response in six degrees of freedom. 
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The accuracy of this process is heavily dependent on the location of the 
a-set. There are analytical techniques available that help find the 'best' a-set 
but ordinarily these are used when applying model reduction methods to find 
modeshapes. These techniques should to be investigated to prove if they can 
also be applied to this process. 

The motivation for this thesis will be realized when the time history 
data of the DDG-51 class shock trial is loaded into I-DEAS, expanded to full- 
order, and then animated in a transient response visualization. 
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APPENDIX A. EXPANDING TIME HISTORIES IN I-DEAS 



Step-by-step procedures for expanding transient time history data using 
I-DEAS software is provided. This option is not available using I-DEAS 
Master Series 3 if the software is set up to run on Hewlett-Packer computers. 
A. ACCESS I-DEAS SIMULATION INTERFACE SOFTWARE 

1. Master Modeler/Surfacing Task - load/build FE model 

2. Meshing Task - mesh model and enter material properties 

3. Boundary Condition Task - 

a. Choose normal mode dynamics 

b. Apply restraints 

c. Create boundary condition set/Choose normal mode 
dynamics - Guyan 

d. Choose analysis set (a-set) 

4. Model Solution Task 

a. Create solution set 

b. Initiate solve 

c. Access I-DEAS model response software 

5. Model Response Software 

a. Excitation definition 

- Create force/displacement excitation function in time 
domain 
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b. Response Evaluation 

- Create response set 

- initiate solve 

c. Save transient response at a-set dof to a file called 'name.unv' 

B. ACCESS I-DEAS TEST INTERFACE SOFTWARE 

1. Correlation Task 

a. import 'name.unv' file and change to an Attached Data File 
(ADF) called 'name, ash' 

b. select mode shape to work with 

- choose 'name.ash' 

- select substructure component 

2. Post-Processing Task 

a. start animation 

note: if using experimental data, once solve is initiated in the Simulation 
Interface Software, the Guyan transformation matrix is created. Proceed 
directly to Test Interface and load the .unv file of a-set time histories 
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APPENDIX B. MATLAB CODE 



A. MAIN FE PROGRAM 

% 

% Program will solve for the system mass and stiffness 

% matrices of a plate. It uses four-node isoparametric 

% elements of the shear-deformable displacement 

% formulation. The system mass and stiffness matrices 

% are saved to a file called sys_mat.mat 
% 

% Written by Scott Waltermire 

% 

% 

% 

% Based on EX 1071. m 

% Prof. Young Kwon 

% 

% 

% 

% 

% 

% 

% Variable descriptions 
% k = element matrix 

% kb = element matrix for bending stiffness 
% ks = element matrix for shear stiffness 
% f = element vector 
% kk = system matrix 
% ff = system vector 
% disp = system nodal displacement vector 
% gcoord = coordinate values of each node 
% nodes = nodal connectivity of each element 

% index = a vector containing system dofs associated with each element 
% pointb = matrix containing sampling points for bending term 
% weightb = matrix containing weighting coefficients for bending term 
% points = matrix containing sampling points for shear term 
% weights = matrix containing weighting coefficients for shear term 
% bcdof = a vector containing dofs associated with boundary conditions 
% bcval = a vector containing boundary condition values associated with 
% the dofs in 'bcdof 

% kinmtpb = matrix for kinematic equation for bending 
% matmtpb = matrix for material property for bending 
% kinmtps = matrix for kinematic equation for shear 
% matmtps = matrix for material property for shear 

% 

% 

% 
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% 

% input data for control parameters 

% 

clear 



hplate=100; 

lplate=100; 

hele=20; 

lele=20; 

a_ele=hele*lele; 

nel=25; 

nnel=4; 

ndof=3; 

nnode=36; 

sdof=nnode*ndof; 

edof=nnel*ndof; 

emodule=30e6; 

rho=0.284/386.4; 

poisson=0.3; 

t=0.125; 

nglxb=2; nglyb=2; 
nglb=nglxb*nglyb; 
nglxs=l; nglys=l; 
ngls=nglxs*nglys; 



% height of plate 
% length of plate 
% element height 
% element length 
% element area 
% number of elements 
% number of nodes per element 
% number of dofs per node 
% total number of nodes in system 
% total system dofs 
% degrees of freedom per element 
% elastic modulus 
% density 
% Poisson's ratio 
% plate thickness 

% 2x2 Gauss-Legendre quadrature for bending 
% number of sampling points per element for bending 
% lxl Gauss-Legendre quadrature for shear 
% number of sampling points per element for shear 



% 

% input data for nodal coordinate values 
% gcoord(ij) where i->node no. and j->x and y 
% size(gcoord) = num_nodes x 2 

% 

% note: if the elements used are not rectangular 
% in shape, the values of gcoord and nodes 

% must be manually entered by the user 

% 



count=0; 
x = 0:lele:lplate; 
y = 0:hele:hplate; 

gnode = zeros(length(y),length(x)); 
gcoord=zeros(length(x) *length(y),2) ; 



for n= 1: length (y) 
for j=l:length(x) 
count = count+l; 

gcoord(count,l:2) = gcoord(count, 1 :2) + [x(j) y(n)]; 
gnode(n,j) = gnode(n,j)+count; 
end 
end 
end 
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% 

% nodes is the nodal connectivity of the model 
% applied in the counter-clockwise direction 

% 

nodes=[]; 

count=0; 

for n=l:length(y)-l 
for j=l .length (x)-l 
count = count+1; 

node(count,l:4)=[gnode(n,j) gnode(n,j+l) gnode(n+l ,j+l) gnode(n+l,j)]; 
nodes = [nodes;node(count,l:4)]; 
end 
end 



% 

% initialization of matrices and vectors 
% 



ff=zeros(sdof,l); % 

kk=zeros(sdof,sdof); % 

mm = zeros(sdof,sdof); % 

disp=zeros(sdof, 1 ) ; % 

index=zeros(edof,l); % 

kinmtpb=zeros(3,edof); % 

matmtpb=zeros(3,3); % 

kinmtps=zeros(2,edof); % 

matmtps=zeros(2,2); % 



system force vector 
system k matrix 
system m matrix 
system displacement vector 
index vector 

kinematic matrix for bending 
constitutive matrix for bending 
kinematic matrix for shear 
constitutive matrix for shear 



% 

% computation of element matrices and vectors and their assembly 

% 

% for bending stiffness 
% 

[pointb,weightb]=feglqd2(nglxb,nglyb); % sampling points & weights 

matmtpb=fematiso(l ,emodule,poisson)*t A 3/12; % bending material property 
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% 

% for shear stiffness 
% 

[points,weights]=feglqd2(nglxs,nglys); 
shearm=0.5*emodule/( 1 .0+poisson); 
shcof=5/6; 

matmtps=shearm*shcof*t*[l 0; 0 1]; 



foriel=l:nel % 

for i=l:nnel 

nd(i)=nodes(iel,i); % 

xcoord(i)=gcoord(nd(i),l); % 

ycoord(i)=gcoord(nd(i),2); % 

end 

k=zeros(edof,edof); % 

m=zeros(edof,edof); % 

kb=zeros(edof,edof); % 

ks=zeros(edof,edof); % 



% sampling points & weights 
% shear modulus 
% shear correction factor 
% shear material property 

loop for the total number of elements 



extract connected node for (iel)-th element 
extract x value of the node 
extract y value of the node 



initialize element stiffness matrix to zero 
initialize element mass matrix to zero 
initialization of bending matrix to zero 
initialization of shear matrix to zero 



% 

% numerical integration for bending term 
% 



for intx=l:nglxb 
x=pointb(intx,l); 
wtx=weightb(intx,l); 
for inty=l:nglyb 
y=pointb(inty,2); 
wty=weightb(inty,2) ; 


% sampling point in x-axis 
% weight in x-axis 

% sampling point in y-axis 
% weight in y-axis 


[shape, dhdr,dhds]=feisoq4(x,y); 


% compute shape functions and 
% derivatives at sampling point 


jacob2=fejacob2(nnel,dhdr,dhds,xcoord,ycoord) ; 

detjacob=det(jacob2); 

in vj acob=inv(j acob2); 


% compute Jacobian 
% determinant of Jacobian 
% inverse of Jacobian matrix 


[dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob); 

% physical coordinate 


% derivatives w.r.t. 


kinmtpb=fekinepb(nnel,dhdx,dhdy); 


% bending kinematic matrix 



% 

% compute bending and mass element matrix 
% 
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kb=kb+kinmtpb' * matmtpb *kinmtpb* wtx * wty *detj acob; 



% 

% Create a matrix of shape functions: Q 

% Ensure each row vector in the matrix is properly weighted: H(i) 

% 

Q = [shape(l) 0 0 shape(2) 0 0 shape(3) 0 0 shape(4) 0 0;... 

0 shape(l) 0 0 shape(2) 0 0 shape(3) 0 0 shape(4) 0; 

0 0 shape(l) 0 0 shape(2) 0 0 shape(3) 0 0 shape(4)]; 

HI = shape(l)*Q; 

HI = [l/12*rho*t A 3*Hl(l,:);l/12*rho*t A 3*Hl(2,:);rho*t*Hl(3,:)]; 
H2 = shape(2)*Q; 

H2 = [l/12*rho*t A 3*H2(l,:);l/12*rho*t A 3*H2(2,:);rho*t*H2(3,:)]; 
H3 = shape(3)*Q; 

H3 = [l/12*rho*t A 3*H3(l,:);l/12*rho*t A 3*H3(2,:);rho*t*H3(3,:)]; 
H4 = shape(4)*Q; 

H4 = [l/12*rho*t A 3*H4(l,:);l/12*rho*t A 3*H4(2,:);rho*t*H4(3,:)]; 

% Create the elemental mass matrix in terms of shape functions 
H = [H1;H2;H3;H4]; 

% 

% Solve for the elemental mass matrix 

% 

m = m+H*wtx*wty*detjacob; 
end 

end % end of numerical integration loop for bending term 

% and mass term 



% 

% numerical integration for shear term 
% 



for intx=l:nglxs 
x=points(intx,l); 
wtx=weights(intx, 1); 
for inty=l:nglys 
y=points(inty,2); 
wty=weights(inty,2) ; 

[shape, dhdr,dhds]=feisoq4(x,y); 



% sampling point in x-axis 
% weight in x-axis 

% sampling point in y-axis 
% weight in y-axis 

% compute shape functions and 
% derivatives at sampling point 



jacob2=fejacob2(nnel,dhdr,dhds,xcoord,ycoord); % compute Jacobian 
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detjacob=det(jacob2); % determinant of Jacobian 

invjacob=inv(jacob2); % inverse of Jacobian matrix 

[dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob); % derivatives w.r.t. 

% physical coordinate 

kinmtps=fekineps(nnel,dhdx,dhdy, shape); % shear kinematic matrix 



% 

% compute shear element matrix 
% 

ks=ks+kinmtps'*matmtps*kinmtps*wtx*wty*detjacob; 

end 

end % end of numerical integration loop for shear term 

% 

% compute element matrix 
% 

k=kb+ks; 



index=feeldof(nd,nnel,ndof); % extract system dofs associated with element 

kk=feasmbl 1 (kk,k,index); % assemble element matrices into system 

mm=feasmbll (mm, m, index); % matrices 

end 



% Apply Springs in z-direction 

kk(3,3) = kk(3,3) + 100; 
kk(18,18) = kk(18,18) + 100; 
kk(93,93) = kk(93,93) + 100; 
kk(108,108) = kk(108,108) + 100; 

save sys_mat mm kk 
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B. FUNCTIONS CALLED BY MAIN FE PROGRAM 



function [point2,weight2]=feglqd2(nglx,ngly) 

% 

% Purpose: 

% determine the integration points and weighting coefficients 
% of Gauss-Legendre quadrature for two-dimensional integration 

% 

% Synopsis: 

% [point2,weight2]=feglqd2(nglx,ngly) 

% 

% Variable Description: 

% nglx - number of integration points in the x-axis 
% ngly - number of integration points in the y-axis 
% point2 - vector containing integration points 
% weight2 - vector containing weighting coefficients 
% 

% determine the largest one between nglx and ngly 

if nglx > ngly 
ngl=nglx; 
else 

ngl=ngly; 

end 



% initialization 

point2=zeros(ngl,2); 

weight2=zeros(ngl,2); 

% find corresponding integration points 



[pointx,weightx]=feglqd 1 (nglx); 
[pointy ,weighty]=feglqd 1 (ngly); 

% quadrature for two-dimension 

for intx=l:nglx 
point2(intx, 1 )=pointx(intx); 
weigh t2(intx, 1 )=weightx(intx) ; 
end 

for inty=l:ngly 
point2(inty,2)=pointy(inty); 
weight2(inty,2)=weighty(inty); 
end 



and weights 

% quadrature rule for x-axis 
% quadrature rule for y-axis 

% quadrature in x-axis 
% quadrature in y-axis 
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function [matmtrx]=fematiso(iopt,elastic, poisson) 



% 

% Purpose: 

% determine the constitutive equation for isotropic material 

% 

% Synopsis: 

% [matmtrx]=fematiso(iopt,elastic,poisson) 

% 

% Variable Description: 

% elastic - elastic modulus 

% poisson - Poisson's ratio 

% iopt=l - plane stress analysis 

% iopt=2 - plane strain analysis 

% iopt=3 - axisymmetric analysis 

% iopt=4 - three dimensional analysis 



if iopt==l % plane stress 
matmtrx= elastic/(l -poisson* poisson)* ... 

[1 poisson 0; ... 
poisson 1 0; ... 

0 0 (l-poisson)/2]; 

elseif iopt=2 % plane strain 
matmtrx= elastic/((l+poisson)*(l-2*poisson))* ... 
[(1 -poisson) poisson 0; 
poisson (1 -poisson) 0; 

0 0 (l-2*poisson)/2]; 

elseif iopt==3 % axisymmetry 
matmtrx= elastic/((l+poisson)*(l-2*poisson))* ... 
[(1 -poisson) poisson poisson 0; 
poisson (1 -poisson) poisson 0; 
poisson poisson (1-poisson) 0; 

0 0 0 (l-2*poisson)/2]; 

else % three-dimension 
matmtrx= elastic/((l+poisson)*(l-2*poisson))* ... 
[(1-poisson) poisson poisson 0 0 0; 

poisson (1-poisson) poisson 0 0 0; 

poisson poisson (1-poisson) 0 0 0; 

0 0 0 (l-2*poisson)/2 0 0; 

0 0 0 0 (l-2*poisson)/2 0; 

0 0 0 0 0 (l-2*poisson)/2]; 

end 
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function [shapeq4,dhdrq4,dhdsq4]=feisoq4(rvalue, svalue) 



% 

% Purpose: 

% compute isoparametric four-node quadilateral shape functions 
% and their derivatves at the selected (integration) point 
% in terms of the natural coordinate 

% 

% Synopsis: 

% [shapeq4,dhdrq4,dhdsq4]=feisoq4(rvalue,svalue) 

% 

% Variable Description: 

% shapeq4 - shape functions for four-node element 
% dhdrq4 - derivatives of the shape functions w.r.t. r 

% dhdsq4 - derivatives of the shape functions w.r.t. s 

% rvalue - r coordinate value of the selected point 

% svalue - s coordinate value of the selected point 

% 

% Notes: 

% 1st node at (-1,-1), 2nd node at (1,-1) 

% 3rd node at (1,1), 4th node at (-1,1) 

% 

% shape functions 

shapeq4( 1 )=0. 25 * ( 1 -rvalue)*( 1 -svalue) ; 
shapeq4(2)=0.25*( l+rvalue)*( 1-svalue); 
shapeq4(3)=0. 2 5 * ( 1 +rvalue) * ( 1 +svalue) ; 
shapeq4(4)=0.25 *( 1 -rvalue) *( 1 +s value); 

% derivatives 

dhdrq4( 1 )=-0.25*( 1 -svalue); 
dhdrq4(2)=0.25*( 1-svalue); 
dhdrq4(3)=0.25 *( 1 -t-svalue) ; 
dhdrq4(4)=-0.25*( 1+svalue); 

dhdsq4(l)=-0.25*(l-rvalue); 
dhdsq4(2)=-0.25*(l+rvalue); 
dhdsq4(3)=0.25 * ( 1 +rvalue) ; 
dhdsq4(4)=0. 25 *( 1 -rvalue) ; 
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functi on [j acob2]=f ej acob2(nnel ,dhdr,dhds ,xcoord,y coord) 

% 

% Purpose: 

% determine the Jacobian for two-dimensional mapping 

% 

% Synopsis: 

% [jacob2]=fejacob2(nnel,dhdr,dhds,xcoord,ycoord) 

% 

% Variable Description: 

% jacob2 - Jacobian for one-dimension 
% nnel - number of nodes per element 

% dhdr - derivative of shape functions w.r.t. natural coordinate r 

% dhds - derivative of shape functions w.r.t. natural coordinate s 

% xcoord - x axis coordinate values of nodes 

% ycoord - y axis coordinate values of nodes 

% 



j acob2=zeros(2,2) ; 
for i=l:nnel 

j acob2( 1 , 1 )=j acob2( 1 , 1 )+dhdr(i) *xcoord(i) ; 
jacob2( 1 ,2)=jacob2( 1 ,2)+dhdr(i)*ycoord(i); 
jacob2(2,l)=jacob2(2,l)+dhds(i)*xcoord(i); 
jacob2(2,2)=jacob2(2,2)+dhds(i)*ycoord(i); 
end 



function [dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob) 

% 

% Purpose: 

% determine derivatives of 2-D isoparametric shape functions with 
% respect to physical coordinate system 

% 

% Synopsis: 

% [dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob) 

% 

% Variable Description: 

% dhdx - derivative of shape function w.r.t. physical coordinate x 

% dhdy - derivative of shape function w.r.t. physical coordinate y 

% nnel - number of nodes per element 

% dhdr - derivative of shape functions w.r.t. natural coordinate r 

% dhds - derivative of shape functions w.r.t. natural coordinate s 

% invjacob - inverse of 2-D Jacobian matrix 
% 



for i=l:nnel 

dhdx(i)=invjacob( 1 , l)*dhdr(i)+invjacob( 1 ,2)*dhds(i); 
dhdy(i)=invjacob(2, 1 )*dhdr(i)+invjacob(2,2)*dhds(i); 
end 
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function [kinmtpb]=fekinepb(nnel,dhdx,dhdy) 



% 

% Purpose: 

% determine the kinematic matrix expression relating bending curvatures 
% to rotations and displacements for shear deformable plate bending 

% 

% Synopsis: 

% [kinmtpb]=fekinepb(nnel,dhdx,dhdy) 

% 

% Variable Description: 

% nnel - number of nodes per element 
% dhdx - derivatives of shape functions with respect to x 

% dhdy - derivatives of shape functions with respect to y 

% 

for i=l:nnel 
il=(i-l)*3+l; 
i2=il+l; 
i3=i2+l; 

kinmtpb( 1 ,i 1 )=dhdx(i); 
kinmtpb(2,i2)=dhdy(i); 
kinmtpb(3 ,i 1 )=dhdy(i) ; 
kinmtpb(3 ,i2)=dhdx(i) ; 
kinmtpb(3,i3)=0; 
end 



function [kk]=feasmbll(kk,k,index) 

% 

% Purpose: 

% Assembly of element matrices into the system matrix 

% 

% Synopsis: 

% [kk] =feasmbl 1 (kk,k, index) 

% 

% Variable Description: 

% kk - system matrix 

% k - element matri 

% index - d.o.f. vector associated with an element 
% 



edof = length(index); 
for i=l:edof 
ii=index(i); 
for j=l:edof 
jj=index(j); 

kk(ii,jj)=kk(ii,jj)+k(i,j); 

end 

end 
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function [kinmtps]=fekineps(nnel,dhdx,dhdy, shape) 

% 

% Purpose: 

% determine the kinematic matrix expression relating shear strains 
% to rotations and displacements for shear deformable plate bending 

% 

% Synopsis: 

% [kinmtps]=fekineps(nnel,dhdx,dhdy,shape) 

% 

% Variable Description: 

% nnel - number of nodes per element 
% dhdx - derivatives of shape functions with respect to x 

% dhdy - derivatives of shape functions with respect to y 

% shape - shape function 

% 

for i=l:nnel 
il=(i-l)*3+l; 
i2=il+l; 
i3=i2+l; 

kinmtps( 1 ,i 1 )=-shape(i) ; 
kinmtps( 1 ,i3)=dhdx(i); 
kinmtps(2,i2)=-shape(i) ; 
kinmtps(2,i3)=dhdy(i); 
end 

function [index]=feeldof(nd,nnel ,ndof) 

% 

% Purpose: 

% Compute system dofs associated with each element 

% 

% Synopsis: 

% [index]=feeldof(nd,nnel,ndof) 

% 

% Variable Description: 

% index - system dof vector associated with element "iel" 

% nd - nodal vector whose system dofs are to be determined 
% nnel - number of nodes per element 
% ndof - number of dofs per node 
% 



edof = nnel*ndof; 
k=0; 

for i=l:nnel 
start = (nd(i)-l) :f: ndof; 
for j=l:ndof 
k=k+l; 

index(k)=start+j; 

end 

end 
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C MAIN POST-PROCESSING PROGRAM 



clear 

% 

% 

% This program will load the system mass and stiffness matrices 
% from file sys_mat.mat. It solves the eigenvalue problem 
% to get the natural frequencies and the modeshapes. 

% Convolution is used to solve for the modal solution q(t). 

% It solves for x(t) at all nodes then calls set.mat to 

% pick out the response at the analysis set. From this it 

% back-expands the aset into the full solution using Guyan and 
% IRS transformation matrices. It can plot the modeshapes, 

% time-histories of each solution and the error between 
% the actual response and the Guyan/IRS models. 

% 

% Written by Scott Waltermire 

% Get the global mass and stiffness matrices 

load sys_mat; 
sdof=108; 

% Create Blast Forcing Function 

plotit = 1 ; 

Fo = 2000; 
time = 0:.01:3; 

[f_of_t] = fBlastForcing(Fo,time,plotit); 

F = zeros(sdof, length (time)); 

F(63,:) = f_of_t; 



% plotit=l plots forcing function 
% amplitude of forcing function 
% time step 

% function to solve for f(t) 

% initialize system force vector 
% place f(t) at the proper node 



% Solve and Sort Eigenvalues/Modeshapes 
[phi,lam]=eig(mm\kk); 

mtilda=phi'*mm*phi; % mass normalize 

% loop to mass normalize the modeshapes 

for i=l :length(mm) 
phi(: ,i)=phi( : ,i) * l/(sqrt(mtilda(i,i))); 
end 

% loop to pick out the natural frequencies 

for j= 1 :length(mm) 
ev(j)=lamO‘,j); 
end 

[lam,p]=sort(ev); % sort the natural frequencies 
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phistar=phi; 

% loop to sort the eigenvectors in the same manner 
% as the natural frequencies 

for k=l;length(mm) 
phi(:,k)=phistar(:,p(k)); 
end 

wn=real(sqrt(lam)) ; 

omega=real((sqrt(lam)/(2*pi)))'; % convert wn to hertz 

phi = real(phi); 

% create modal_F and use proportional damping 

modal_F = phi'*F; 
zeta = .02; 

% Convolution to solve for q(t) using Trapezoidal rule 

dt= .001; 
for i= 1:108 

wd(i) = wn(i)*sqrt(l-zeta A 2); 
h = l/wd(i)*exp(-zeta*wn(i)*time).*sin(wd(i)*time); 

F = modal_F(i,;); 

[convXY] = fconvTrap(F,h,dt); 
q(i,:) = convXY; 
end 

% Plot Modeshapes 

num_mode_shape= 1 0 ; 
for i=l:num_mode_shape 
phi_mode=real(phi(:,i)); 
phi 1 =(phi_mode(3 : 3:18))’; 
phi2=(phi_mode(2 1 :3:36))’; 
phi3=(phi_mode(39:3:54))'; 
phi4=(phi_mode(57:3:72))’; 
phi5=(phi_mode(7 5:3:90))’; 
phi6=(phi_mode(93 : 3 : 1 08))' ; 
phi_mesh=[phi 1 ;phi2 ;phi3 ;phi4 ; phi5 ;phi6] ; 
surfc(phi_mesh) 
shading interp 
grid 
pause 
end 
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% Convert back to physical Coordinates 



resp = phi*q; 



% Solve for T_static and T_irs 

load set; % aset/oset file 

% functions to solve for Guyan and IRS 

% transformation matrices. 

[kstat,mstat,T_static] = fstatic_tam(kk,mm,oset,aset); 
[kirs,mirs,T_irs] = firs_tam(kk,mm,oset,aset); 



% 

% Load Guyan Transformation matrix 
% and determine xa and xo 
% 

for i = lilength(aset) 
xa_g(i,:) = resp(aset(i),:); 
end 

% Back expand the guyan aset time histories 
% into the full response 

x_g = T_static*xa_g; 

xo_g = x_g(l:sdof-length(aset),:); 

xa_g = x_g(sdof-length(aset)+l:sdof,:); 

% Now re-sort xa_g and xo_g 

for i=l :length(aset) 
xol_g = xo_g(l:aset(i)-l,:); 
xo2_g = xo_g(aset(i):length(xo_g(:,l )),:); 
xol_g = [xol_g;xa_g(i,;)]; 
xo_g = [xol_g;xo2_g]; 
end 

x g = xo_g; 



%• 



% Load IRS Transformation matrix 
% and determine xa and xo 



%■ 



for i = 1 :length(aset) 
xa_irs(i,:) = resp(aset(i),:); 
end 
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% Back expand the irs aset time histories 
% into the full response 

x_irs = T_irs*xa_irs; 

xo_irs = x_irs(l:sdof-length(aset),:); 

xa_irs = x_irs(sdof-length(aset)+l:sdof,:); 

% Now re-sort xa_irs and xo_irs 

for i= 1 :length(aset) 
xol_irs = xo_irs(l:aset(i)-l,:); 
xo2_irs = xo_irs(aset(i):length(xo_irs(:,l)),:); 
xol_irs = [xol_irs;xa_irs(i,:)]; 
xo_irs = [xol_irs;xo2_irs]; 
end 

x_irs = xo_irs; 



% Ignore all rotations and plot response 
% in the z-direction 

resp = resp(3:3:sdof,:); 
resp_g = x_g(3:3:sdof,:); 
resp_irs = x_irs(3:3:sdof,:); 

for i=l:36 

figure(l) 

plot(time,resp(i,:),'r-.’,time,resp_g(i,:),y); 
xlabel('Time (sec)'); 
ylabel(’Displacement (in)'); 
title('Time History'); 

figure(2) 

plot(time,resp(i,:)/r-.',time,resp_irs(i,:),'y'); 
xlabel('Time (sec)'); 
ylabel('Displacement (in)'); 
title('Time History'); 

pause 

end 

close(l) 

close(2) 
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% 

% Find the error between the true response and 
% the Guyan/IRS back-expanded response 
% 

for i = l:length(time) 

xg_err(:,i) = abs(resp_g(:,i) - resp(:,i)); 
xirs_err(:,i) = abs(resp_irs(:,i) - resp(:,i)); 
end 

fori= 1:36 

xg_error(i) = max(xg_err(i,:)); 
xirs_error(i) = max(xirs_err(i,:)); 
end 

xg_errorl = xg_error(l:6); 
xg_error2 = xg_error(7:12); 
xg_error3 = xg_error(13:18); 
xg_error4 = xg_error( 19:24); 
xg_error5 = xg_error(25:30); 
xg_error6 = xg_error(31:36); 

xg_error = [xg_errorl;xg_error2;xg_error3;xg_error4;xg_error5;xg_error6]; 

figure(l) 

surf(xg_error); 

shading interp 

grid 

xirs_error 1 = xirs_error(l:6); 

xirs_error2 = xirs_error(7:12); 

xirs_error3 = xirs_error(13:18); 

xirs_error4 = xirs_error( 19:24); 

xirs_error5 = xirs_error(25:30); 

xirs_error6 = xirs_error(31:36); 

xirs_error = [xirs_errorl;xirs_error2;xirs_error3;... 

xirs_error4;xirs_error5;xirs_error6]; 

figure(2) 
surf(xirs_error); 
shading interp 
grid 



save modal_info resp resp_g resp_irs time 
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D. FUNCTIONS CALLED BY MAIN PROGRAM 

function [f_of_t] = fBlastForcing(Fo,time,plotit); 

% 

% This function returns a forcing function which is 
% a "blast" function. 

% 

% F(t) = Fo * ( exp(-at) - exp(-bt) ) 

% 

% where a and b are constants which shape the blast, 

% and Fo is the amplitude of the blast. 

% 

% The variable "plotit" is a switch which if set = 1 will 
% cause the f(t) to be plotted, if set to anything else 
% will not plot. 

% written by J.H. Gordis 

% 

% Choices: sine blst step 

type = 'blst'; 

% 

if type == 'blst'; 



disp(' Blast forcing used...') 
a = 100.0; 
b = 300.0; 

f_of_t = Fo * ( exp(-a*time) - exp(-b*time) ); 



elseif type == 'step'; 



% disp(' Step forcing used...’) 
f_of_t = Fo * ones(size(time)); 



elseif type == 'sine'; 



% disp(' Sine forcing used...’) 

W = 5; % Hertz 

f_of_t = Fo * sin(2*pi*W*time); 



end; 



if plotit == 1 ; 
figure(9) 

plot(time,f_of_t);grid 

pause 

elf 

end 

% End function. 
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function [convXY] = fConvTrap(x,y,dt); 

% 

% Usage: [convXY] = fConvTrap(x,y,dt); 

% 

% This function performs the convolution integral calculation 

% 

% tau=t 

% convXY(t) = S x(t-tau) y(tau) d tau 
% tau=0 

% 

% using the trapezoidal rule. 

% The function is passed the vectors x(t) and y(t) 

% and the sample spacing (time step) dt. 

% 

% x = vector of length n 

% y = vector of length n 

% dt = sample spacing (i.e. delta_t). 

% written by J.H. Gordis 

% 

if size(x) == size(y); % Vectors the same size. Perform convolution. 

% 

convXY = zeros(size(x),l); 

for icnt_step = 2 : length(x); 

[wts] = fTrapzWts(icnt_step); 
if size(x,l) -= size(wts,l); 

wts = wts'; 
end 

convXY(icnt_step) = dt * sum(wts .* fliplr(x(l:icnt_step)) .* y(l:icnt_step) ); 
end; 

elseif size(x) ~= size(y) & length(x) = length(y); % Vectors are transposed. Don't 
perform convolution. 

% 



disp(' Error in fConvTrap. Transpose one vector and try again.') 

disp(sprintf(' Size(x) %3i',size(x))) 
disp(sprintf(' Size(y) %3i',size(y))) 

elseif size(x) ~= size(y) & length(x) -= length(y); % Vectors different sizes. Don't 
perform convolution. 

% 



disp(' Error in fConvTrap. Vector not the same size.') 

disp(sprintf(' Size(x) %3i',size(x))) 
disp(sprintf(' Size(y) %3i',size(y))) 



end; 
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% 

function [kstat,mstat,T_static]=fstatic_tam(k,m,oset,aset) 

% 

% this function returns the statically reduced stiffness 
% and mass matrices, given the unreduced couterparts. 

% Care must be taken that the aset and oset vectors correspond 
% with the existing arrangement of k and m. 

% k and m are UNPARTITIONED matrices. 

% 

% written by J.H. Gordis 

aset_size=length(aset); 

% 

kaa=k(aset,aset); 
kao=k(aset,oset); 
koo=k(oset,oset); 
koa=kao'; 
clear k; 

k=[koo,koa;kao,kaa]; 

% 

maa=m(aset,aset); 
mao=m(aset,oset); 
moo=m(oset,oset); 
moa=mao'; 
clear m; 

m=[moo,moa;mao,maa] ; 

% 

t_static=-koo\koa; 

T_static = [t_static;eye(aset_size)]; 

% 

kstat=T_static'*k*T_static; 

mstat=T_static' :t: m :i: T_static; 

% 

% end function fstatic tarn 
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% 

function [kirs,mirs,T_irs]=firs_tam(k,m,oset,aset) 

% 

% this function returns the IRS reduced stiffness 
% and mass matrices, given the unreduced couterparts. 

% Care must be taken that the aset and oset vectors correspond 
% with the existing arrangement of k and m. 

% k and m are UNPARTITIONED matrices. 

% written by J.H. Gordis 

aset_size=length(aset); 

% 

kaa=k(aset,aset); 
kao=k(aset,oset); 
koo=k(oset,oset); 
koa=kao’; 
clear k; 

k=[koo,koa;kao,kaa]; 

% 

maa=m(aset,aset) ; 
mao=m(aset,oset); 
moo=m(oset,oset); 
moa=mao'; 
clear m; 

m= [moo,moa;mao ,maa] ; 

% 

t_static=-koo\koa; 

T_static = [t_static;eye(aset_size)]; 

% 

kstat=T_static'*k*T_static; 

mstat=T_static'*m*T_static; 

% 

tirs=t_static+inv(koo)*(moa+moo*t_static)*inv(mstat)*kstat; 

T_irs= [tirs;eye(aset_size)] ; 

% 

kirs=T_irs'*k*T_irs; 

mirs=T_irs'*m*T_irs; 

% 

% end function firs tarn 
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File holding the dof for the analysis set (aset) 
and the omitted set (oset) 

Saved to a file called set. mat 



% 

% 

% 

% 

% 

% 

% 

%■ 



aset = [3 18 93 108]; 

oset = [1 2 4 5 6 7 8 9 10 1 1 12 13 14 15 16 17 ... 

19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 ... 

40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 ... 
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 ... 

83 84 85 86 87 88 89 90 91 92 94 95 96 97 98 99 100 101 102 ... 

103 104 105 106 107]; 

save set aset oset 



% Program will animate the full-order response 
% by loading a .mat file called modal-info 

clear; 

load modaHnfo; 
clear resp_g resp_irs 

% make movie 

movie_fig = figure('position’,[100 200 300 200]); 
M = moviein(200); 

[x,y] = meshgrid([-5:2:5]); 
for i= 1:200 
resp_t = resp(:,i); 
respl = (resp_t(l:6))'; 
resp2 = (resp_t(7:12))'; 
resp3 = (resp_t(13:18))’; 
resp4 = (resp_t( 19:24))'; 
resp5 = (resp_t(25:30))'; 
resp6 = (resp_t(31:36))'; 
z = [respl ;resp2;resp3;resp4;resp5;resp6]; 
surf(x,y,z); 
shading interp 
grid; 

axis([-5 5 -5 5 -.4 .6]); 
view([45 45 10]) 
zlabel('disp (in.)'); 

M(:,i) = getframe; 

end 

pause 

movie(M,0); 
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% Program will animate the plate response found 
% by the static reduction method. It loads the .mat 
% file called modal_info 

clear; 

load modal_info; 

clear resp resp_irs 
resp = resp_g; 

% make movie of guyan response 

movie_fig = figure('position',[100 200 300 200]); 

M = moviein(200); 

[x,y] = meshgrid([-5:2:5]); 

count = 1 ; 
for i= 1:200 
resp_t = resp(:,i); 
respl = (resp_t( 1:6))'; 
resp2 = (resp_t(7:12))'; 
resp3 = (resp_t(13:18))'; 
resp4 = (resp_t( 19:24))'; 
resp5 = (resp_t(25:30))'; 
resp6 = (resp_t(31:36))'; 
z = [respl;resp2;resp3;resp4;resp5;resp6]; 
surf(x,y,z); 
shading interp 
grid; 

axis([-5 5 -5 5 -.4 .6]); 
view([45 45 10]) 
zlabel('disp (in.)'); 

title(’ Animation by Guyan Back Expansion'); 

M(:, count) = getframe; 
count = count+1 ; 
end 

pause 

movie(M,0); 
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% 

% Program to animate the plate response found 

% by the IRS reduction method. Loads the .mat file 
% modal_info 

% 

clear; 

load modal_info; 

clear resp resp_g 
resp = resp_irs; 

% make movie of guyan response 

movie_fig = figure('position',[100 200 300 200]); 

M = moviein(200); 

[x,y] = meshgrid([-5:2:5]); 

count = 1 ; 
for i = 1 ;200 
resp_t = resp(:,i); 
respl = (resp_t(l:6))'; 
resp2 = (resp_t(7:12))'; 
resp3 = (resp_t(13:18))'; 
resp4 = (resp_t(19:24))'; 
resp5 = (resp_t(25:30))'; 
resp6 = (resp_t(31;36))'; 
z = [respl;resp2;resp3;resp4;resp5;resp6]; 
surf(x,y,z); 
shading interp 
grid; 

axis([-5 5 -5 5 -.4 .6]); 
view([45 45 10]) 
zlabel('disp (in.)’); 

title('Animation by IRS Back Expansion'); 

M(:, count) = getframe; 
count = count+1; 
end 

pause 

movie(M,0); 
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